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ABSTRACT 

We study the distribution function (DF) of dark matter particles in haloes of mass 
range 1O^''-1O^^M0. In the numerical part of this work we measure the DF for a sample 
of relaxed haloes formed in the simulation of a standard ACDM model. The DF is 
expressed as a function of energy E and the absolute value of the angular momentum 
L, a form suitable for comparison with theoretical models. By proper scaling we obtain 
the results that do not depend on the virial mass of the haloes. We demonstrate that 
the DF can be separated into energy and angular momentum components and propose 
a phenomenological model of the DF in the form fE{E)[l + L2/(2L2)]-/3oo+/3ox-2/3o_ 
This formulation involves three parameters describing the anisotropy profile in terms 
of its asymptotic values (/3o and Poo) and the scale of transition between them (Lo)- 
The energy part fsiE) is obtained via inversion of the integral for spatial density. 
We provide a straightforward numerical scheme for this procedure as well as a simple 
analytical approximation for a typical halo formed in the simulation. The DF model is 
extensively compared with the simulations: using the model parameters obtained from 
fitting the anisotropy profile, we recover the DF from the simulation as well as the 
profiles of the dispersion and kurtosis of radial and tangential velocities. Finally, we 
show that our DF model reproduces the power-law behaviour of phase space density 
Q^p{r)/a^{r). 
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1 INTRODUCTION 

The distribution function (DF) provides the most general 
and complete way of statistical description of dark matter 
(DM) haloes. It carries maximum information on the spatial 
and velocity distributions of particles in such objects. Our 
knowledge on the DF is still being improved, mostly due to 
numerical experiments. In the last few years cosmological 
simulations have revealed increasingly detailed features of 
phase-space structure of DM haloes. These numerical results 
provide useful constraints on theoretical models of the DF. 
One property of interest in this field is the anisotropy of the 
velocity dispersion tensor. It has been demonstrated that the 
outer parts of the haloes exhibit more radially anisotropic 
trajectories than the halo centre (see e.g. Coh'n, Klypin & 
Kravtsov 2000; Fukushige & Makino 2001; Wojtak et al. 



2005; Mamon & Lokas 2005; Cuesta et al. 2007). This fea- 
ture, besides the well-studied density profile, has been con- 
sidered as the main point of reference in the attempts at 
construction of a reliable model of the DF. 

So far, a few approaches to this problem have been pro- 
posed. Cuddeford (1991) generalized the Osipkov-Merritt 
model (Osipkov 1979; Merritt 1985) to the DF which gener- 
ates an arbitrary anisotropy in the halo centre and becomes 
fully radial at infinity. Although an analytical inversion for 
these models exists, the anisotropy profile cannot be rec- 
onciled with the numerical results: the rise from central to 
outer anisotropy is too sharp and the outer orbits are too 
radial (see Mamon & Lokas 2005). An & Evans (2006a) no- 
ticed that a non-trivial profile of the anisotropy can be ob- 
tained from a sum of DFs with a constant anisotropy for 
which an analytical inversion is known (Cuddeford 1991; 
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Kochanek 1996; Wilkinson & Evans 1999). However, the re- 
sulting anisotropy profiles axe decreasing functions of ra^ 
dius and do not agree with those measured in cosmologi- 
cal simulations. Recently a very elegant method has been 
presented by Baes & van Hese (2007). The authors intro- 
duced a general ansatz for the anisotropy profile and then, 
for a given potential-density pair, derived the DF as a series 
of some special functions. This approach works well under 
the condition that the potential can be expressed as an ele- 
mentary function of the corresponding density. This require- 
ment, however, is not satisfied by many models, including 
the NFW density profile (Navarro, Frenk & White 1997) 
which is commonly used as a good approximation of the 
universal density profile of DM haloes. 

The DF inferred from the simulation gives a possibil- 
ity to test directly the analytical models. According to the 
Jeans theorem any spherically symmetric system in the state 
of equilibrium should possess a DF which is a function only 
of energy and the absolute value of angular momentum. This 
theoretical postulate was taken into account in the compu- 
tation carried out by Voglis (1994) and Natarajan, Hjorth & 
van Kampcri (1997) . In the first case the DF was obtained for 
a single relaxed halo which formed from cosmologically con- 
sistent initial conditions. It was shown that there were two 
main contributions to the DF, the halo population and the 
core population of particles. Both were effectively described 
by two independent phcnomenological fits. Natarajan et al. 
(1997) determined the DF for a sample of cluster-size haloes 
formed in cosmological simulations. Their selection of ob- 
jects included those with substructures and departing from 
equilibrium. They also discussed and took into account in 
their calculation the effect of boundary conditions defined by 
the virial sphere. However, the final results were not used to 
test quantitatively any model of the DF. 

It seems that two main approaches to study the DF, 
namely theoretical modelling and feedback from the simula- 
tions, evolved rather separately barely crossing each other. 
The rare exceptions include the work of Lokas & Mamon 
(2001) who used the Eddington formula to derive numeri- 
cally the DF following from the NFW profile in the isotropic 
case and that of Widrow (2000) who considered more general 
cuspy profiles and Osipkov-Merritt anisotropy. This paper is 
devoted to combining both approaches and providing a co- 
herent analysis of the DF from the viewpoints of the simula- 
tions as well as the model construction. Our main aim is to 
propose a phcnomenological model of the DF that recovers 
the results from the simulations as accurately as possible. 

Our effort is mainly motivated by the future appli- 
cations of the derived DF to the dynamical modelling of 
galaxy clusters. Although subhaloes in general have differ- 
ent density and velocity distributions than DM particles 
(Diemand, Moore & Stadel 2004), massive subhaloes (those 
likely to host galaxies) are distributed like DM particles (Fal- 
tenbacher & Diemand 2006). Although the correspondence 
between the massive subhaloes and galaxies in real clusters 
remains to be proven, our results should at least in principle 
be applicable to kinematic data sets for galaxy clusters. The 
traditional approach to do such modelling was to reproduce 
the velocity dispersion profile of the galaxies by solving the 
Jeans equation (see e.g. Katgert, Biviano & Mazure 2004). 
It is well known however that from the dispersion alone 
one cannot constrain all the interesting parameters (such 



as the virial mass, the concentration of the NFW profile 
and anisotropy) because of the density-anisotropy degener- 
acy. One can break this degeneracy by using the fourth or- 
der velocity moment, the kurtosis, and solving an additional 
higher-order Jeans equation (Lokas 2002; Lokas & Mamon 
2003; Lokas et al. 2006; Wojtak & Lokas 2007). Although 
this approach has many advantages (e.g. it does not require 
the knowledge of the full DF) , it has been applied till now 
ordy for constant-anisotropy models and the calculation of 
velocity moments requires the binning of the data in which 
some information is lost. Since the number of galaxies with 
measured redshifts per cluster is still rather low (of the order 
of a few hundred for the best-studied, nearby clusters) it is 
essential that all the available information is used. This can 
be obtained by fitting the projected DF to the data directly. 

A few approaches along these lines have been attempted 
already. For example, Mahdavi & Geller (2004) used a simple 
DF of the form f{E,L) oc ^"-1/2^-2/3 (^^ich yields con- 
stant anisotropy) to constrain the mass profile and orbital 
structure using combined kinematic data sets for nearby 
galaxy groups and clusters, while van der Marel et al. (2000) 
in their study of CNOC clusters did not assume an explicit 
form for the energy-dependent part of the DF, but still used 
the constant anisotropy. Wojtak et al. (2007) used a sim- 
plified form of the projected isotropic DF constructed from 
the projected density combined with a Gaussian distribu- 
tion for the linc-of-sight velocities to study the properties 
of members versus interlopers in simulated kinematic data 
sets. None of the DFs used so far, however, reflects accu- 
rately the true properties of cluster-size DM haloes found in 
TV-body simulations. 

The paper is organized as follows. Section 2 provides 
the theoretical framework and defines all the basic quanti- 
ties used later on in the paper. In the next section we discuss 
the details of the computation of the DF of DM particles in 
the haloes formed from cosmological simulations and provide 
examples of the results. Section 4 is devoted to the deriva- 
tion of a phcnomenological model of the DF; we discuss the 
separability of the DF in energy and angular momentum 
and present an explicit formula for the L-dcpendent part of 
DF. An extensive comparison of the model with the simu- 
lations is presented in Section 5, where we also provide an 
analytical approximation for the energy-dependent part of 
the DF obtained for an average halo. Finally, the discussion 
follows in Section 6. 



2 THEORETICAL FRAMEWORK 

This section summarizes the theoretical background of the 

paper. First, we introduce scaling properties consistent with 
the NFW density profile. We will use this profile in the pa- 
per, but our approach is not restricted to this particular 
density distribution and can be easily generalized to any 
profile consistent with simulations (see below). Second, we 
briefly describe the relation between the differential DF and 
the DF itself. Finally, we discuss the consequences of the fi- 
nite volume of the virialized area of a halo. In particular, the 
relation between the DF and its differential form is properly 
modified to account for this effect. 
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2.1 Scaling properties 

It is a well known fact that the density profiles of DM haloes 
formed in cosmological simulations exhibit striking similar- 
ity. NFW showed that most of them are well fitted within the 
virial sphere by the universal two-parameter profile which 
can be expressed in the following way 

, . 1 M. 1 

where x — r /r^. The two free parameters are the scale radius 
ra and the mass enclosed within the sphere of this radius 
Ms. The (positive) gravitational potential inferred from the 
Poisson equation reads (Cole & Lacey 1996; see also Lokas 
& Mamon 2001) 



vE'NPw(r) = K 



;ln(l -f x) 



(2) 



where the velocity unit Vs is related to the circular velocity 
V^ir(''a) at the scale radius via Vs = Vcir(f's)(ln2 — 1/2)""'^''^. 

Let us note that rs, Vs and Ms define a set of natu- 
ral units of the NFW model. By scaling any quantity by 
a proper combination of them we remove the explicit de- 
pendence on the free parameters of the NFW model. This 
is an essential property if we want to study the dynamics 
of a whole class of haloes with NFW-like density profiles. 
Hereafter, we will keep this scaling in all equations in the 
text. In many places we will also use a unit of the angular 
momentum Ls as a substitute for V^rg. 



2.2 The distribution function 

The DF is a fundamental concept in statistical mechanics 
of A''-body systems. It describes the phase-space density of 
particles of such a system without any detailed knowledge 
of the time evolution of A*' trajectories. Following the Jeans 
theorem, a steady-state DF, which is of interest for us here, 
depends on the phase-space coordinates only through the 
integrals of motion. Although the shape of DM haloes is in 
general better approximated by a three-axial ellipsoid rather 
than a sphere (see e.g. Gottlober & Yepes 2007), it is still 
very efi'ective to assume spherical symmetry in dynamical 
approach. Given that the streaming motions and internal 
rotation within the virial sphere are negligible compared to 
higher velocity moments, spherical symmetry implies that 
the DF can be expressed as 



/(r,v)=/(iJ,L), 



(3) 



where E is the positively defined binding energy and L the 
absolute value of the angular momentum per unit mass 



E = *(r) - iv^ 
^ ' 2 

L = r X V. 



(4) 
(5) 



The gravitational potential in equation Q is related to the 
DF through the Poisson equation 



(6) 



V^*(r) = -47rG / /(r,v)d''v 



The most natural and straightforward probe of f{E, L) 
in numerical experiments is the so-called diff'erential DF de- 
fined in the following way 



N{E,L) = 



(7) 



dEdL 

One may intuitively interpret this function as mass density 
in energy-angular momentum space. The DF itself can be 
simply derived dividing N{E,L) by the volume g{E,L) of 
the hypersurface of constant energy and angular momentum 
embedded in the phase space 

N{E,L) 



f[E,L) = 



9{E,L) 



(8) 



It is easy to show that the volume of this hypersurface reads 
(see Appendix A) 



g{E,L) = 8n' LTriE,L), 



(9) 



where Tr{E,L) is the radial period of an orbit given by 
the following integral over radius from the pericentre to the 
apocentre 

dr 



Tr 



(10) 



rp V2[*('-) - -B - iV(2r2)] ■ 

The upper panel of Fig. [1] shows a contour map of 
g{E, L) (dotted lines) calculated for the NFW gravitational 
potential ([2]). The I/max(-E) line is the profile of maximum 
angular momentum which consists of points corresponding 
to circular orbits. This curve divides the energy-angular mo- 
mentum plane into an area describing the physical orbits of 
a system (below Lmax) and the zone not permitted by me- 
chanics (above Lmax). Note that we are using the scaling 
relations introduced in the previous subsection so the re- 
sults do not depend explicitly on the halo mass Ms and the 
scale radius r^. In some places later on we will refer to the 
inverse function for Lmax(i5) by i5max(i). 

Voglis (1994) and Natarajan et al. (1997) showed that 
the dependence of Tr{E,L) on the angular momentum is 
very weak and could be neglected without loss of precision. 
This is understandable if we note that the NFW-like po- 
tentials are still not so far away from the isochrone poten- 
tial ^'(r) oc (6 -I- \/b'^ + r'^)~^ which leads to purely energy- 
dependent Tr proportional to E~^^^ (Binney & Tremaine 
1987). Following Natarajan et al. (1997) we will use this 
feature to simplify expression To do this we first note 
that the volume of the hypersurface of constant energy qe 
is given by 

gEiE) = r'^^^^^\iE,L)dL. (11) 
Jo 

Taking advantage of the weak dependence of Tr on L, equa- 
tion (|10|) . we get 



g{E,L)^2 



9e{E)L 



(12) 



-^'max(-£^) 

On the other hand, one can show that gsiE) reads (see 
Appendix A) 



gE{E) = 167r^ / V2(*(r) - E^dr, 



(13) 



where rma.^{E) is the apocentre radius of the radial orbit. In- 
serting (|13|l into (|12|) one immediately gets a very simple ap- 
proximation for g{E, L) involving only a one-dimensional in- 
tegral without singularities, in contrast with g{E, L) derived 
by expression (|10|l . We find that this approximation repro- 
duces the exact formula (O with enough accuracy. Taking 
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Figure 1. The volume of the hypersurface of constant energy 
and angular momentum for the infinite halo (upper panel) and 
the finite halo limited by the virial sphere of radius = Sr^ 
(lower panel). In both cases the NFW profile was assumed. Solid 
lines show the profiles of the maximum angular momentum of a 
given system. In the lower panel the three shades of gray mark 
the three characteristic zones according to the orbit size defined 
by the relation of the virial radius Tv to the pericentre radius rp 
and the apocentre radius ra, as labelled. 



advantage of its numerical simplicity we use it in majority 
of our calculations. 



2.3 Boundaries of the haloes 

So far we have discussed the relation between the DF and its 
differential form for an infinite system. In practice, however, 
we restrict our numerical analysis to the interior of the virial 
sphere which separates the equilibrium part of a halo from 
the infall region. We define the virial radius of this sphere 
by 



4 3, 

-TrrvAcPc 



(14) 



where Mv is the virial mass, pc is the present critical density 
and Ac is the virial overdensity. Another parameter com- 
monly used to describe the size of the virial sphere in terms 
of rs is the concentration c = r^/vs. 

The existence of the boundary of the virialized part of 
the halo implies that Tr given by (|10|) must be replaced by 



, , (15) 

V2[*(r)-S-L7(2r2)] 

where the upper limit of the integral is a minimum of the 
virial radius and the radius at the apocentre (see Appendix 
A for details). Combining p5|l with one gets a general 
formula for the volume g{E, L) in the presence of a spherical 
boundary of a halo. Contrary to the conclusion of Natara- 
jan et al. (1997), we find that the approximation (I12|) is 
no longer justified for orbits extending beyond the virial 
sphere (ra > r^). This follows from the fact that angular 
momentum dependence of (|15p becomes non-negligible and 
the integral (|11|) cannot be simplified to the form of (|12|l . 

Using ((9)1 and (fTS)) with the NFW potential, we cal- 
culated g{E,L) for a halo limited by the virial sphere of 
radius = br^ (see the lower panel of Fig. [T|. As expected, 
the result differs from an infinite system by the orbits with 
Ta > rv and remains unchanged for trajectories wholly in- 
cluded within the virial sphere. 



3 THE DISTRIBUTION FUNCTION FROM 
THE SIMULATION 

3.1 The simulation 

For our TV-body simulation we have assumed the WMAP3 
cosmology (Spergel et al. 2007) with matter density f2,n — 
0.24, the cosmological constant SIa = 0.76, the dimension- 
less Hubble parameter h = 0.73, the spectral index of pri- 
mordial density perturbations n = 0.96 and the normaliza- 
tion of power spectrum erg = 0.76. We have used a box of 
size 160/i~'^Mpc and 1024^ particles. Thus the particle mass 
was 3.5 X 10* Mq . Starting at redshift 2 = 30 we followed the 
evolution using the MPI version of the Adaptive Refinement 
Tree (ART) code (Kravtsov, Klypin & Khokhlov 1997). 

We identified clusters with the hierarchical friends-of- 
friends (FOF) algorithm (Klypin et al. 1999) with a linking 
length of 0.17 times the mean inter-particle distance which 
roughly corresponds to an overdensity of 330. We have se- 
lected 36 clusters at redshift 2 = in the range of virial 
mass (0.15-2) x IO^^Mq, where the virial overdensity param- 
eter appropriate for our cosmological model was assumed 
Ac = 93.8 (Lokas & Hoffman 2001). Our sample did not in- 
clude clusters with two substructures of approximately the 
same mass and a poor fit of the NFW profile suggestive of 
a recent major merger. 

Starting from the FOF position of the cluster we have 
determined the highest density peak as the final centre of the 
clusters. This centre coincides with the position of the most 
massive substructure found at the linking length 8 times 
shorter and also with the position of the halo found by the 
EDM halo finder (Klypin et al. 1999). 

3.2 Computation of the distribution function 

In the first step of the computation we calculate the binding 
energy (|4| and the angular momentum ([5} per unit mass for 
each particle within the virial sphere of each halo. Spherical 
symmetry implies that we have to apply in (2} the radial 
profile of the gravitational potential 

GM{r)dr 



^(r) 



f 



(16) 
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Figure 2. Energies and angular momenta of particles witliin the 
virial sphere of one of the simulated haloes. Solid and dashed 
lines mark the profile of the maximum angular momentum and 
the limit of vanishing radial velocity at rv respectively. To make 
the picture less obscure we plotted only 1 percent of the particles. 



where >I'(oo) = 0. However, the mass profile of the equihb- 
rium part of a halo reaches no further than the virial radius. 
On the other hand, all analytical models of the DF involve 
the density profile extending to infinity. We found that the 
only coherent way to reconcile both facts is to split the in- 
tegral (|16l) in two parts 



*(r) = 



GM{r)dr 



+ 



GMKFw(r)dr 



(17) 



The first term is evaluated numerically by integration of a 
discrete mass profile. The second term is an analytical exten- 
sion with the NFW density profile which is an assumption of 
the DF model introduced in the following section. Its contri- 
bution to the potential is a constant equal to ln(l + c) /c. 

Fig. [2] shows the resulting energies and angular mo- 
menta of particles inside the virial sphere of one of the sim- 
ulated haloes. The profile of the maximum angular momen- 
tum (solid line) and the profile of vanishing radial velocity 
at the virial sphere (dashed line) were calculated for the 
exact gravitational potential given by (|17p . All particles oc- 
cupy the area permitted by mechanics or lie very close to the 
boundary line. Interestingly, quite a large fraction of them 
have orbits extending beyond the virial sphere. As noted in 
the previous section, we keep and Lb as units of energy 
and angular momentum respectively. The parameters of the 
NFW model were obtained for each halo by fitting the NFW 
formula to the density profile measured in logarithmic radial 
bins. 

In the next step we determine for each halo the dif- 
ferential DF given by ([Hi. In this calculation we used our 
own version of the FiEstAS (Field Estimator for Arbitrary 
Spaces) algorithm designed to infer the density field from a 
scatter diagram embedded in a space of any number of di- 
mensions (see Ascasibar & Binney 2005 for more details). As 
a result of this computation we get an estimate of N{E, L) 
at all points of the energy- angular momentum plane cor- 
responding to the particles inside the virial sphere. Once 
N[E, L) is calculated the DF can be easily obtained via (|8]). 
As discussed in section 2, we used approximation (|12p for 
the orbits contained inside the virial sphere and the exact 
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Figure 3. Contour maps of the DF of DM particles inside the 
virial sphere of two example haloes. The profile of the maximum 
angular momentum is indicated by Lmax(^) and the line of van- 
ishing radial velocity at the virial sphere by ra = r-v (apocentre 
at the virial sphere) or rp = (pericentre at the virial sphere). 



formula ((9)l with (|15p for trajectories extending beyond r\. 
We found that the additional advantage of expression (|12|) is 
that it could be evaluated at any point of the energy-angular 
momentum plane. This helps us to keep the estimates of 
the DF obtained for points with angular momentum lying 
slightly above Lmax(S). 

In order to derive a contour map or a profile of the DF 
we introduce a regular dense mesh on the energy-angular 
momentum plane and find the median value of the DF in 
each cell. Such a set of median points is considered as the fi- 
nal numerical approximation of the DF and is used in prepa- 
ration of all plots in this paper. Fig. [3] shows two exam- 
ples of the resulting contour maps obtained for two different 
haloes. The unit of the DF in this and following Figures is 
Ms/r^ /Vs ■ The interval between the iso-DF lines is fixed 
at value 0.25 of the logarithmic scale. The lack of the DF 
estimation in the lower part of each diagram arises from the 
fact that this zone is occupied by very few particles (see e.g. 
Fig. [5)) so that no information on the distribution can be 
retrieved. Let us note that this is an effect of the finite mass 
resolution. 
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Figure 4. The profile of the anisotropy parameter. The sofid hne 
and the dark gray area are the median and the interquartile range 
of the profiles obtained for individual haloes and rescaled by 
inferred from fitting the NFW profile. 



4 THE ANALYTICAL MODEL OF THE 
DISTRIBUTION FUNCTION 

A general form of the DF for spherical systems in the state 
of equilibrium is a function of energy and the absolute value 
of angular momentum f{E, L). In our approach we assume 
that the DF is separable in energy and angular momentum 



f{E,L) = fEiE)fLiL). 



(18) 



This is the first assumption that considerably narrows the 
family of possible solutions. Therefore, it is necessary to 
check how robust it is. We address this problem in the next 
section, where we present an extensive comparison of the 
analytical model with the simulations. 

The angular momentum part of the DF in equation 
UlSp specifies the anisotropy of velocity dispersion tensor. 
This quantity is commonly described with the so-called 
anisotropy parameter 



1 



(19) 



where ct^ and ag are the radial and the tangential velocity 
dispersions respectively and we assume there are no stream- 
ing motions. The values of this parameter range from — oo 
for circular orbits to 1 for purely radial trajectories. Fig. 2] 
shows the average anisotropy profile of the simulated haloes 
used for the measurement of the DF. The fight gray rect- 
angle in the background of the plot indicates the position 
of the virial radius. It is clearly seen that the anisotropy is 
typically a growing function of radius, with values ~ 0.07 
in the halo centre and ~ 0.3 at the virial sphere (see e.g. 
Mamon & Lokas 2005 and Cuesta et al. 2007 for compari- 
son). On the other hand, the considerable width of the in- 
terquartile range of the measured /3(r) (dark gray region) 
signifies that the profiles of single haloes differ among each 
other. Occasionally flat or decreasing profiles are measured. 
It seems that a simple and general enough analytical model 
of the anisotropy should possess at least three free parame- 
ters which determine asymptotic values of P{r) for small and 
large radii and a scale of transition between them. We pro- 
ceed with the construction of such a model by introducing 
a proper ansatz for /z, (L) . 



Louis (1993) showed that the following asymptotes of 
the angular momentum part of the DF 



for L < io 
for L > Lo, 



(20) 



where Lo is an angular momentum constant, lead to con- 
stant anisotropy /3oo at infinity (r^*I'(r) 3> Lq) and (3 = 
in the halo centre. This result can be easily generalized to 
the case of a non-isotropic velocity distribution in both lim- 
its of radius. First, let us note that the central part of the 
halo is dominated by the particles with small angular mo- 
menta, namely < 2r^^{r) <^ Lq. Then, remembering 
that the DF of constant anisotropy takes the form (Henon 
1973; Binney & Tremaine 1987; Lokas 2002) 

f{E,L)=fE{E)L-^^, (21) 

it is easy to notice that the formula (|20|l can be rewritten in 
the following way 

/ L-^^° forL«Lo 

where l3o is the central anisotropy of a system. As shown by 
An & Evans (2006b), the upper limit for f3o is equal to 7/2, 
where is the density profile near the halo centre. This 
means that for the NFW density model we have /3o < 1/2. 

The simplest function obeying the asymptotic condi- 
tions formulated above is a double power-law function 

-/3oo-l-/3o 



-2/3o 



(23) 



As shown in the following section, this ansatz leads to a 
very realistic anisotropy profile that fits well the l3{r) profiles 
of simulated haloes. Furthermore, the simplicity of formula 
P3|) guarantees that the energy part of the DF can be quite 
easily calculated via the inversion of the integral equation 

The key idea of this procedure lies in an analytical simpli- 
fication of the right-hand side of (|24p to a one-dimensional 
integral. The resulting equation is then solved numerically 
for fsiE). The technical details of this calculation are sum- 
marized in Appendix B. Once the full form of the DF is 
determined one can also calculate the velocity moments. All 
formulae are reduced to one-dimensional integrals which can 
be easily evaluated numerically (see Appendix C). 

The top row of Fig. [5] shows the anisotropy, disper- 
sion (Tr and kurtosis Kr = {v^)/a^ of the radial velocity 
inferred from the model of the DF. The calculations were 
carried out assuming the NFW density profile and four sets 
of model parameters chosen to illustrate the flexibility of the 
model: /3o — 0.1 and /3oo ~ 0.3, 0.5 (solid and dashed lines 
respectively); /3o = /9oo = 0.3 (dotted line); /3o = 0.4 and 
Poa = 0.1 (dashed-dotted line). In all cases the transition 
value of Lo = 0.25 Ls was used. 

The dispersion profiles for the two models with increas- 
ing /3(r), as expected, differ only for large radii which is the 
effect of different values of Paa ■ Interestingly, the correspond- 
ing kurtosis profiles clearly signify fiat-topped velocity distri- 
bution in the outer part of the halo (k^ < 3), highly peaked 
distribution in the centre (k,. > 3) and roughly Gaussian for 
radii around {nr ~ 3). On the other hand, non- increasing 
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Figure 5. The top panels show (from left to right) the anisotropy, dispersion and kurtosis of the radial velocity inferred from the model 
of the DF with four sets of parameters: /3o = 0.1, Lo = 0.25 La and jiaa = 0.3, 0.5 (solid and dashed lines respectively); /3o = /3tx3 = 0.3, 
Lq = 0.25 Ls (dotted line); /3o = 0.4, Lq = 0.25 Ls and /3oo = 0.1 (dashed-dotted line). The corresponding DFs for the same sets of 
parameters are plotted in the bottom panels in terms of: the energy part of the DF fsiE) (left), iso-DF lines with values indicated 
along the curve of maximum angular momentum (middle) and the profiles of the DF for three values of angular momentum (right). In 
all calculations the NFW density profile was assumed. 



/3(r) profiles lead to less peaked velocity distributions in the 
centre. It seems therefore that the typical anisotropy of DM 
haloes, as shown in Fig. |4l is expected to coincide with the 
kurtosis rapidly growing towards the halo centre (see also 
Fig. [TO] below). As we will see in the following section, this 
is one of the most characteristic features of the phase-space 
structure of massive DM haloes. 

In the bottom panels of Fig. [S]we plotted the DFs cor- 
responding to four sets of model parameters, as described 
above. The three panels from the left to the right show the 
energy part of the DF, contour maps and the profiles for 
three fixed values of angular momentum. The plots reveal 
some interesting signatures of the specific shape of /3(r) pro- 
file. For example, the inclination of the iso-DF lines with re- 
spect to the energy axis decreases with increasing f3o'- more 
isotropic /3 at the centre corresponds to more vertical iso- 
DF lines; also the shape of the lines is somewhat different. 
These features are also to some extent visible in the contour 
maps of the DF for two simulated haloes in Fig. |3] The up- 
per map represents a halo with an increasing f3{r), whereas 
the second one depicts the case of a decreasing f3{r) profile. 
Both haloes are analyzed in terms of velocity moments and 
the DF in the following section. 

Recently Mamon & Lokas (2005) found that the simu- 
lation data are well reproduced by the anisotropy profile of 
the form 




Figure 6. Comparison between the functional form of the 
anisotropy I I25I I proposed by Mamon & Lokas (2005) (/3ml) 
and /3(r) inferred from the DF model (/3op) with /3o = and 
Poc = 0.5. 



Fig. [6] shows both anisotropies for three values of ri/4. Note 
that both j3{r) profiles have similar shapes, although our 
anisotropy profile has a somewhat sharper rise at small radii. 

We also find that the radius ro characteristic of the DF 
model, for which /3 is the mean of the limiting values 



/3M = 



2r + ri/4' 



(25) 



/?o + /3o 



(26) 



where ri/4 is the radius where (3 = 0.25. Assuming l3o — 
and /3oo = 0.5 in our DF model, we made a comparison 
of the resulting anisotropy with the functional form (|25p . 



depends weakly on Po and jSoo ■ For parameter ranges leading 
to P{r) profiles covering the interquartile area of anisotropy 
from the simulation (0 < /3o < 0.15, < /3oo < 0.6 and 
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0.04 < L(i/Lb < 25), this radius is well (within 5 percent 
accuracy) approximated by 

ro/r, = 3.69(Lo/Ls)°-^^ + 2.27[Lo/L,f-\ (27) 



5 COMPARISON WITH THE SIMULATION 

5.1 The distribution function 

The DF proposed in the previous section is a phenomeno- 
logical model in the sense that it possesses free parameters 
whose values should be adjusted to the simulation data. All 
three parameters were introduced to determine a family of 
anisotropy profiles so that it is /3(r) that is most sensitive 
to the variations of /3o, Poo and Lq. Consequently, we de- 
cided to constrain the parameters of the model by fitting 
the P{r) profile inferred from the DF model to the median 
profile measured in simulated DM haloes. The best-fitting 
parameters are: /3o = 0.09, /3oc = 0.34 and Lq = 0.198 Lg. 
The corresponding best-fitting profile of the anisotropy is 
plotted as a dashed line in the lower left panel of Fig. [9] 

Once the model parameters are adjusted the DF can be 
compared with its counterpart measured from the simula- 
tion. Fig.[7]shows this comparison in terms of a contour map 
and the profiles for constant angular momentum or energy. 
Dark gray regions in all panels indicate the interquartile 
range of the DF values within the halo sample. The lighter 
gray area in the background of the upper diagram marks the 
points of vanishing radial velocity at the virial radius . Its 
boundaries are fixed by the first and third quartile of virial 
radii in the halo sample, 4.1ra and 6. Org respectively. 

Although some deviations of the model (dashed lines) 
from the results of the simulations are visible, in general 
the theoretical profiles are included within the interquartile 
range or lie very close to its boundaries. As expected, the 
strongest discrepancy between the model and the simula- 
tion is present in the part of the energy- angular momentum 
plane populated by the particles with orbits extending be- 
yond the virial sphere (the area to the left of the r^, = 
line). However, given that this is the only part of the energy- 
angular momentum plane affected by the infalling material, 
we think that the observed differences are acceptable. 

5.2 The separability of the distribution function 

A critical point of the derivation of the DF presented in the 
previous section was the factorization introduced by equa- 
tion (|18p . In order to inspect the robustness of this assump- 
tion we propose a simple test. We calculate the ratio of the 
DF from the simulation to the energy part of the DF model 
with parameters adjusted to the anisotropy profile from the 
simulation. Under the assumption that the real DF is factor- 
izable in energy and angular momentum, we can expect that 
the resulting ratio should be a weak function of energy equal 
to fL(L) given by (|23|) . Fig. [8] shows that the variations of 
this ratio with respect to fL{L) are of the same order as 
the width of the interquartile range which means that sep- 
arability is acceptable from the statistical point of view. A 
small systematic deviation can be seen for L ~ 0.1 Ls. How- 
ever, this is certainly a local feature since this trend is not 
repeated in other profiles. Let us emphasize that this test 
of separability depends strongly on the reliability of /_l(L). 
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Figure 7. Comparison between the DF measured for DM parti- 
cles inside the virial sphere of the simulated haloes and the model 
with parameters adjusted to the median anisotropy, /3o = 0.09, 
Poo = 0.34 and Lq = 0.198 Ls. Solid lines and gray areas stand 
for median profiles and interquartile ranges of the DF measured 
in the halo sample, whereas the dashed lines correspond to the 
model. The light gray area in the background of the upper dia- 
gram indicates the points of vanishing radial velocity at the virial 
radius Tv. Its boundaries are fixed by the first and third quartile 
of virial radii in the halo sample. 



One can imagine that an incorrect form of /z, (L) would likely 
lead to a negative result of the test, whether f{E, L) is sepa- 
rable or not. On the contrary, a positive result of such a test 
in our case means that not only is the assumption of factor- 
ization valid but the approximation for /l(L) is reasonable 
as well. 
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Figure 8. The ratio of tlie DF inferred from the 36 simulated 
haloes to the energy part of the DF model with parameters ad- 
justed to the median anisotropy profile. Each panel shows the 
profile for the constant value of angular momentum indicated in 
the lower left corner. Solid lines and gray areas represent the 
median profiles and interquartile ranges respectively. The dashed 
lines indicate the values of given by the ansatz II23I I. 



5.3 Velocity moments 

Further comparison between the simulation and the DF 
model can be done in terms of velocity moments. This is 
depicted in Fig. [9] where the dispersion and kurtosis of the 
radial and tangential velocity are plotted. In the bottom part 
of this figure we show the profiles of the anisotropy /3(r) and 
/34 parameter which measures the anisotropy of a tensor of 
the fourth velocity moment. By analogy with the parameter 
/3(r) we defined /34(r) in the following way 



Mr) 



(28) 



The dashed lines in each panel of Fig. |9] are the model pre- 
dictions, except for the f3{r) profile (lower left panel) which 
is a fit of the model providing constraints on parameter val- 
ues given in the previous subsection. Theoretical dispersion 
profiles coincide very well with the profiles from the simu- 
lation. We notice quite a good agreement also for the P4(r) 
parameter. On the other hand, theoretical profiles of the 
kurtosis are systematically biased towards higher values, but 
typically by less than 10 percent. Nevertheless, their shapes 
clearly recover the shapes of the median profiles from the 
simulation. Moreover, for both the radial and tangential ve- 
locity a characteristic growth of k from value < 3 around 
the virial radius up to > 4 in the halo centre, is seen. 

Although the kurtosis bias is enclosed within acceptable 
limits (kurtosis is known to be sensitive to any noise), it 
would be desirable to find out the reason for this behaviour. 
Since our statistical samples consist of 10* — 10^ particles 
per radial bin, we ruled out a possibility of a bias of the 
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Figure 9. Velocity moments of the radial and tangential velocity 
(top and middle panels) and the anisotropy of the dispersion ten- 
sor (S{r) and the fourth velocity moment tensor /34(r') defined by 
II28I I. The solid line and the dark gray area are the median and the 
interquartile range of the profiles obtained for individual haloes 
and rescaled by rs and Vb inferred from fitting the NFW profile. 
Dashed lines are the profiles inferred from the DF model with pa- 
rameters adjusted to the median anisotropy measured from the 
simulation. In the case of anisotropy (S{r) (lower left panel) this 
line shows the fit of the model. 



kurtosis estimator (see Lokas & Mamon 2003). We also ex- 
cluded the possibility that this is caused by some specific 
assumptions of the model. For example, changing the NFW 
density distribution to the 3D Sersic profile, which fits the 
simulation data even better (Navarro et al. 2004; Merritt et 
al. 2005; Prada et al. 2006), we still encounter the same bias. 
In addition, perturbing the model parameters of /l(I/) does 
not explain the situation either. 

We therefore conclude that the slight discrepancy in 
the predictions of our model concerning the kurtosis must 
signify reaching the limitations of the theoretical approach 
based on using the global, smooth gravitational potential of 
a system. We suppose that the problem is caused by the 
presence of substructures which perturb locally the trajec- 
tories of particles with respect to the orbits determined by 
the global potential of a halo. What one gets from the sim- 
ulation is really a convolution of the velocity distribution 
expected from the model involving a global potential with 
the distribution of velocity perturbations occurring due to 
density fluctuations. The estimation of the importance of 
this effect is a complicated task since the perturbation of 
the particle orbit depends on many variables, such as the 
distribution of substructures, softening of the potential and 
particle velocity. However, some qualitative conclusions can 
be drawn. First, note that low- velocity particles are affected 
by the density perturbations more strongly. Consequently, 
the peak of the resulting velocity distribution is suppressed 
and the tails are preserved which may effectively decrease 
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the kurtosis (see Fig. [9]). Second, the effect of the perturba- 
tion on the velocity dispersion is a higher order correction 
compared to the dispersion obtained for a system with a 
global potential. This means that the resulting dispersion 
profiles are barely changed and they are still expected to 
coincide well with theoretical predictions. 

It seems an intriguing issue that the profile of tangential 
kurtosis signifies Gaussianity of the velocity distribution at 
radii around where the logarithmic slope of the density 
profile is equal to —2. One could suppose that some signa- 
tures of the so-called isothermal sphere are locally present. 
Interestingly, this statement is also supported by the shape 
of the DF for E < *(rs) 0.7 that is the energy range 
of particles at rs- Referring to the middle panel of Fig. [7] 
it is easy to notice that the DF grows exponentially with 
energy, as expected for systems not very different from the 
isothermal sphere. The distribution of the radial velocity, 
on the other hand, takes the Gaussian form for radii around 
0.3 rs. This difference could be a consequence of the non- 
vanishing anisotropy parameter, which is not accounted for 
in the classical formulation of the isothermal model (Binney 
1982; Binney & Tremaine 1987). 

So far we tested the DF model for a typical halo as- 
sociated with the median properties of our halo sample. In 
order to check the applicability of our model more exten- 
sively we repeat such comparison for single haloes. The DF 
in this case is expected to differ from one halo to another 
due to the observed variety of anisotropy profiles. Results of 
this analysis are summarized in Fig. 1101 To save space we in- 
cluded only five haloes with representative, rather different 
anisotropy profiles (upper panels), from the most strongly 
increasing profile in the left panel to a decreasing one on 
the right-hand side. The second and fifth panels correspond 
to the haloes for which contour maps of the DF are shown 
in Fig. [3] (the top and bottom panel respectively) . We re- 
stricted the number of profiles to those most essential: we 
plot the dispersion and kurtosis of the radial velocity and the 
anisotropics /3(r) and /34(r). We also show the profiles of the 
DF for three values of angular momentum or energy. In all 
panels the solid lines represent simulation results, whereas 
dashed lines are the predictions of the model. As in Fig. |9l 
dashed lines in the case of parameter /3(r) indicate best fit- 
ting profiles of the model. Gray regions in the panels of two 
bottom rows mark the interquartile ranges of the DF which 
describe the scatter of points resulting from the FiEstAS 
algorithm. 

From the analysis of Fig. [TU] we conclude that all pro- 
files, regardless of the anisotropy, are very well reproduced 
by our model of the DF. In general, the theoretical DF does 
not exceed the limits of the interquartile range or lies very 
close to its boundaries (see two bottom rows of panels in 
Fig. llOf) . Surprisingly, we find that the agreement is usu- 
ally almost equally good when the model is applied to the 
haloes with massive substructures which were rejected from 
our sample. This is certainly good news for the future ap- 
plications of our DF to the dynamical modelling of galaxy 
clusters which very often display signatures of recent major 
mergers. 




Figure 11. Relative errors of velocity moments and anisotropies 
inferred from the DF obtained with the analytical approximation 
of (E) given by i'29t . All profiles were compared with the results 
of exact calculations summarized in Appendix C. 



parameter 


/3o = 0.09 
/3oo = 0.34 


/3o = 
l3ao = 0.5 


C 


0.000738 


0.001663 


Eo 


0.12903 


0.14124 


El 


0.078 


0.07 


ai 


2.10 


1.74 


E2 


0.071 


0.085 


02 


2.47 


3.01 



Table 1. Values of the parameters used in the approximation 
of the energy part of the DF 12911 . The first column lists the 
parameters. The second column gives the parameter values for 
the model fitted to the anisotropy for the average halo from the 
simulation (/3o = 0.09 and /3oo = 0.34) and the third one gives 
the values which reproduce the DF for the anisotropy profile II25II 
(/3o = and /3oo = 0.5). In both cases Lq = 0.198 Ls was assumed. 



5.4 Analytical approximation of the distribution 
function 

The DF discussed in the first subsection is typical in a sense 
that it describes the statistical macrostate of DM particles 
in a typical massive halo. With future applications in mind 
we decided to provide an analytical approximation for the 
energy part /b(-E) of this DF which could be used as a 
substitute for a rather complicated procedure described in 
Appendix B. We found that the following expression repro- 
duces the numerical DF with good accuracy 

fE{E) = Cexp('-|-')i5°i/[^+^''P(^''^i)l 



X(l _ ^)-"2/{l+cxp[(l-B)/B2)]}^ 



(29) 



where the values of the parameters are listed in the mid- 
dle column of Table [T] For completeness we recall that the 
angular momentum part of the DF is given by ((23} with 
/3o = 0.09, /3oo = 0.34 and Lo = 0.198 Ls. We have ver- 
ified that the errors of the dispersion, kurtosis and both 
anisotropies, when using this approximate formula in the in- 
tegrals for velocity moments, do not exceed 5 percent within 
the radial range (O.Olr-s, 30rs) (see Fig. [TTJ. Note that the 
general form of expression (|29|l can be effectively used to ap- 
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Figure 10. Comparison between the model of the DF and the simulation for 5 haloes (in columns) in terms of the anisotropics /3(r) and 
/34(r), the dispersion and kurtosis of the radial velocity and the profiles of the DF for different values of angular momentum or energy. 
Solid and dashed lines show respectively the results from the simulation and the predictions of the model with parameters adjusted to 
the anisotropy profile. Gray areas and lines in the panels of two bottom rows indicate the interquartile ranges and the medians of the 
DF for three fixed values of angular momentum or energy given in the upper left corner of each panel (with lower profiles corresponding 
to higher angular momentum or lower energy). 



proximate the DF model also for other sets of parameters. 
As a second example we include in the third column of the 
Table also the parameters of a model with Po — 0, /3oo ~ 0.5 
and I/O = 0.198 Ls which mimics the anisotropy profile (|25p 
with ri/4 = 0.9 rg. 



6 DISCUSSION 

We have studied the DF of DM particles inside the virial 
spheres of the haloes of mass 10"-1015Mq formed in the 
standard ACDM cosmological A^-body simulation. In the 
first part of the paper we presented results of the calcu- 
lation of the DF from the simulation in the form most suit- 
able for comparison with theoretical models. Then we pro- 
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posed a phenomenological model of the DF. The model in 
its part dependent on angular momentum involves three free 
parameters which specify the anisotropy profile, namely its 
asymptotic values and the scale of transition between them. 
We demonstrated that this parametrization is sufficient to 
reproduce accurately the simulation results in terms of ve- 
locity moments as well as the DF itself. The only discrepant 
point we encountered was a small but statistically significant 
bias of the theoretical kurtosis with respect to its profiles 
measured from the simulation. This is probably caused by 
the presence of substructures perturbing the trajectories of 
low-velocity particles. 

In section 5 we showed that the velocity distribution of a 
typical halo changes from a flat-topped distribution (k < 3) 
in the outer part to a peaked one (k > 3) near the centre. 
This behaviour was noticed and discussed by others before 
(e.g. Kazantzidis, Magorrian & Moore 2004; Wojtak et al. 
2005). The analysis of the DF presented here suggests that 
this property of the velocity distribution is correlated with 
the profile of the anisotropy increasing with r: I3{r) profiles 
growing faster with r imply more rapid growth of the kur- 
tosis towards the centre. 

As demonstrated by Taylor & Navarro (2001), the pro- 
file of the phase-space density Q{r) = p{r)/a{r)^ in DM 
haloes is well fitted by a power-law function. It seems that 
the status of this relation is as well established as the NFW 
fit of the density profile. We checked that Q{r) profiles in- 
ferred from the DF model with parameters adjusted to the 
median anisotropy (see section 5) coincide well with the cor- 
responding power- law functions (see Fig. I12p . In this com- 
parison we assumed the logarithmic slopes of —1.92 (in the 
case of the dispersion of the radial velocity) and —1.84 (in 
the case of the total dispersion), the values obtained from 
the simulations by Dehnen & McLaughlin (2005). The rela- 
tive residuals of both Q{r) profiles are of the same order as 
the scatter of points from the simulations in fig. 1 of Dehnen 
& McLaughlin (2005). Note that this happens when the DF 
model is tuned to the mean trend of the P{r) parameter. 
Therefore one could suspect that both relations, the mean 
profile of the anisotropy and Q{r) oc r are two aspects 
of some deeper relation. A more general parametrization of 
the DF might provide some insights towards a more funda- 
mental understanding of this phenomenon. 

Although the whole analysis presented in this paper was 
done in the framework of the NFW density profile, the equa- 
tions for the numerical inversion (IB12|I were derived for an 
arbitrary density distribution. Using this general form one 
can immediately obtain a family of DFs with our general 
anisotropy profiles for any potential-density pair. For the 
commonly used density profiles it is easy to introduce the 
phase space units analogous to rs, 14 and Afs in our case. 
This reduces the role of the parameters of the density profile 
to scaling properties so that the final DF model would not 
explicitly depend on them. 

Given our very general parametrization of the /3(r) pro- 
file, our DF model is expected to provide some impact on the 
solution of the classical problem of mass-anisotropy degen- 
eracy for spherical systems. In order to obtain more reliable 
estimates of mass profiles, one could assume the anisotropy 
profile from the simulation and keep the density profile as 
the only degree of freedom of the DF. One could then apply 
the maximum likelihood approach of projected DF, as de- 




10"- 10"' lo" lo' 10- 10"^ 10"' 10° lo' 10^ 



Figure 12. Radial profiles of p(r)/ar{r)'^ (left panel) and 
p{r) /(Jtotir)^ (right panel). The dashed lines show predictions of 
the model DF with the NFW density profile and parameters fitted 
to the median anisotropy profile shown in Fig. |4] The solid lines 
plot power-law functions with logarithmic slopes from Dehnen & 
McLaughlin (2005). In the lower panels we show relative residuals. 

scribed e.g. by Mahdavi & Geller (2004). A more advanced 
and simulation-independent approach would be to treat the 
anisotropy profile as an unknown quantity, described by the 
three parameters introduced in our formulation. As a result 
one would obtain an estimate of the mass profile as well as 
anisotropy. Both methods require an additional study of the 
DF in projection and extensive tests on mock data sets. This 
will be the subject of our follow-up papers. 
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APPENDIX A: THE VOLUME OF THE HYPERSURFACE 

The volume g{E, L) of the hypersurface Se l(v, r) of constant energy and angular momentum is defined by the integral 
g{E,L)AEAL= ( A^vA^r. (Al) 

isEt(v,r)dBdL 

Introducing spherical coordinates and changing variables into E, L and radius r one gets 

g(E, L)AEAL = Stv'^LAEAL <£ - — — --, (A2) 

where Vr is the radial velocity and the integral is equal to the radial period of the orbit. Using the radii of the pericentre rp 
and the apocentre ra, one can rewrite the final formula for g{E,L) in the following way 

g(E,L) = 16n^L n '^'^ (A3) 

Integrating g{E,L) over the angular momentum we get the volume gE{E) of the hypersurface of constant energy 
gE{E) = 1671^ / LAL / ^ (A4) 

Changing the order of the integrals and performing the integral over the angular momentum one obtains 



gE(E) = 167r^ / \/2(*(r) - E^Ar, (A5) 



where rmax(S) is the apocentre radius of the radial orbit {E — ^'(rmax))- 

Considering a system of finite size V(v, r) in phase space, one has to recalculate g{E, L) with the realistic hypersurface of 
constant E and L given by the Cartesian product Se l(v, r) x V{v, r). In particular, for a spherical system with a boundary 
in the form of a sphere of radius the upper limit of the integral in (|A3|I and (IA5|) must be replaced by min{ra,r-v} and 
min{ra,rmax(-E)} respectively. 



APPENDIX B: THE ENERGY PART OF THE DISTRIBUTION FUNCTION 

The energy part of the DF /_b {E) introduced in section 4 is related to the density profile by 



'"^'^ ^ III ^''^^X^ 2lf j L-'^^>A\. (Bl) 

Although the main part of this paper concerns the DF consistent with the NFW profile, we keep within the appendix a general 
density p(r) so that the final formulae of inversion could be applied to any potential-density pair. 
Changing variables in the integral (|B1[I into the energy and angular momentum one gets 

^(^) = a^/^-^Vr-^Lj-^"" /* fE{E)AE r (^ + ^);!l!!°^" -dA, (B2) 

Jo Jo ^/x — X 

where x = r'^{-^ - E) / Ll and A = L'y{2Ll). The integral over the A variable is evaluated analytically so that (|B2|) can be 
rewritten in the form 

2;3o _ /n^N3/2„-/3o rC!_-A)) r .pw.-r, RNl/2-;io; 



p(r)r^^" = (2^)^^^2-^» / fEm'^~EY''''"K[^,E)AE (B3) 

1 (3/2 - Ih) Jo 

with a kernel of the integral given by 

K{^,E) = (l + a;)-''°°+^«2-Fi(l/2,/3oo -/3o,3/2-/3o,x/(l + a;)), (B4) 



where 2F1 stands for the hypergeometric function. Equation (|B3} is a Volterra integral of the first kind. In the general case of 
models with varying anisotropy, when /3ao 7^ A) and Lo < 00, it has no analytical solution for Je^E) due to the complexity of 
expression (|B4[I . However, as shown by Cuddeford & Louis (1995), this kind of integral can be quite easily inverted numerically. 
Below we adapt their method to our problem. 

For E ^ and 'if ^ E the integral kernel can encounter a singularity, i.e. K(if,E) oc E^°°~^° . In order to avoid this 
behaviour, we define a smooth integral kernel K{'^ , E) which is free of such a feature 

i?(*, E) = £;"''~+'^«it'(*, E). (B5) 
By analogy we introduce a smooth energy part of the DF which is a regular function for energy approaching 
fE{E)=E'''^-''fE{E). (B6) 



The value of v is determined from the limit of p(r) at small potential (^' <C 1 and r 00) in equation (|B3|) . for which 
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p{ry^« oc iSj(ry'^Or^(-f>--+0o)^ (B7) 

which immediately gives v in terms of /3o, /3oo and the asymptotic slopes of the potential and density. In particular, for the 
NFW density profile one gets i/ = 3 — 2Poo + Po- 

Using formulae HB5|I and HB6|) we can rewrite equation (IB3|I in the following form 



Jo 

where p = p{r)r'^^° and Cp^ stands for all coefficients in front of the integral (|B3|l . Following Cuddeford & Louis (1995) we 
introduce discrete vectors of the potential = je, radius rj = r{'^ j) and density pj = p{rj)r'^'^'', where j is an integer number 
and e — 1/j <^ 1. For any pj we can split the integral HB8|I into a sum 

P. = Cp, £ r fE{E)K{je, E){je - £) V2-/3o^-3/2+,3^-,3„^^^ ^gg^ 

In order to apply any numerical algorithm to invert equation (|B9[) with respect to fE{E), one has to assume that e is 
sufficiently small so that the variations of fsiE) and K{je,E) within subsequent integration ranges are negligible. Then one 
can approximate both functions by their values at {i — l/2)e, i.e. the middle points of the energy range. This approach was 
used by Cuddeford & Louis (1995) and favoured over other methods involving higher order interpolation (see e.g. Saha 1992). 
Applying this approximation to equation (|B9P we get 



P. = C0, J2 f^^^U^' - l/2)6)/>, (BIO) 

i = l 

with fsi ~ fE{{i — l/2)e) and the matrix hj defined in the following way 

hj = (ie)''+''~-'''''[B,/.('^ + /3oo-A)-l/2,3/2-/3o)-B(._i)/,(i^ + /3co-A)-l/2,3/2-/3o)], (BU) 



where Bz{x,y) is the incomplete beta function. As shown by Cuddeford & Louis (1995), the solution of (|B10[l for fsi can be 
obtained by evaluating iteratively the following expression 

f _ Pj/Cpo - E-Ij-^ Kjje, ii - l/2)e)IijfEi 
with the initial value /e i given by 

fEi- J'^^'° . (B13) 
if(e,e/2)/n 



APPENDIX C: VELOCITY MOMENTS 

All non- vanishing moments of the radial or tangential velocity at a given radius r are defined by the integral 

1 r r r / r^x -i^tx+Po 

^"^^^^'^ ^ w) IIL'^H' ^ ^-^'-^^"^'^^ (CD 

where i is a subscript indicating the velocity component in spherical coordinates and n is an integer number. By analogy 
with steps (IB2P and (|B3[) from Appendix B, one can evaluate the angular momentum part of the integral and simplify the 
expression to the following form 

93/2+n-/3o /■* 

= \,;3„^(^) fE{E){^ - Ef'+-^°K.{^,E)dE, (C2) 

where the kernel functions for the radial (i = r) and tangential {i = d — (f)) component of the velocity field are given 
respectively by 

Kr{<l',E) = /?y\rri/?i~ ^ + ^)-^-+'''2fi(l/2 + n,p^- /3o, 3/2 + n- /3o, x/(l + x)) (C3) 
(1/2 + n)l (3/2 + n — flo) 

Kei^,E) = If/l + ""^^1 + - /;°) (1 + :c)-^°°+^",Fi(l/2, f3^ - p,, 3/2 + n - p,,x/[l + x)). (C4) 
r(l + n)r(3/2 + n — Po) 



Once /b(-B) is calculated in terms of the fEi vector it is easy to evaluate numerically the integral (|C2|I and obtain the profile 
of any velocity moment. 
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An interesting property of the model is the ratio of any non-vanishing moment of the tangential velocity to the corre- 
sponding moment of the radial velocity in the limit of small and large radii. Introducing spherical coordinates in (|C1|) and 
performing the integral with two asymptotes of /l(-£/) given by (|22p. one can show that this ratio is the following function of 

/3o or Poa 



= < 



r(l-/3o)r(l + n) 
r(l + n-/3oo) 

I. r(i - /3oo)r(i n) 



for r ^ 



(C5) 



for r 
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